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Abstract. We have measured the muon flux at the underground Gran Sasso National 
Laboratory (3800 mw.e.) to be (3.41 ±0.01) • 10 _4 m _2 s _1 using four years of Borexino data. 
A modulation of this signal is observed with a period of (366±3) days and a relative amplitude 
of (1.29 ± 0.07)%. The measured phase is (179 ± 6) days, corresponding to a maximum on 
the 28 th of June. Using the most complete atmospheric data models available, muon rate 
fluctuations are shown to be positively correlated with atmospheric temperature, with an 
effective coefficient ay = 0.93 ± 0.04. This result represents the most precise study of the 
muon flux modulation for this site and is in good agreement with expectations. 
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1 Introduction 

Muons observed in underground sites arise mostly from the decay of pions and kaons produced 
by the interaction of primary cosmic rays with the nuclei of the upper atmosphere [1]. Since 
muons lose energy as they penetrate the Earth, low-energy muons are stopped and only the 
most energetic muons can be observed in underground detectors, with an energy threshold, 
S t hr, depending on the depth. The flux of cosmic muons detected deep underground shows 
time variations which are, at first approximation, seasonal. These variations can be related to 
the air density fluctuations, which affect the fraction of mesons decaying to muons energetic 
enough to reach the underground detector. This effect has been known and studied for 
many decades [2]. Underground experiments have observed this phenomenon at Gran Sasso 
(MACRO [3], LVD [4]) and in other underground sites ([5], [6] and refs. therein). 

Borexino is an organic liquid scintillator detector [7] located in the underground Gran 
Sasso National Laboratory (LNGS) in central Italy under a limestone coverage of ~1300m 
(~3800mw.e). It is devoted to the spectroscopy of low-energy solar neutrinos via elastic 
scattering on electrons. Data taking started in May 2007 and led to measurements of solar 
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neutrinos ( 7 Be [8, 9], pep [10], 8 B [11], and a limit on CNO [10]), as well as antineutrinos 
from the Earth (geo- neutrinos) [12]. Borexino is also a powerful tool for both the study of 
cosmic muons that penetrate the Gran Sasso rock coverage and the neutrons and radioactive 
isotopes that they produce, which are relevant backgrounds for dark matter and neutrino 
experiments. 

Borexino can select muons passing through a spherical volume with a cross section of 
146 m 2 . Such a geometry makes the acceptance independent of the muon angle of incidence, 
allowing us to measure the muon flux and its time dependence with reduced systematics. 
Furthermore, as air temperature data can be obtained from specialized atmosphere modeling 
centers for weather forecast [13], the correlation with the muon flux can be investigated and 
the effective temperature coefficient can be determined. 

In this article we first introduce the Borexino detector (section 2) and report on the 
measured muon flux (section 3) and on its seasonal modulation (section 4). We then briefly 
describe how the muon flux is expected to be related to the atmospheric effective temperature 
(section 5), before reporting the temperature fluctuations at LNGS (section 6). Both muon 
flux and temperature modulations are also analyzed with Lomb-Scargle frequency analy- 
sis (section 7). Finally we report the correlation between the muon flux and atmospheric 
temperature (section 8) before summarizing our results (section 9). 

Preliminary results of this analysis have been presented in [14]. 

2 The Borexino Detector 

The Borexino detector [7] is schematically shown in Figure 1. The active target for the 
analyses reported in this article is the Inner Detector (ID). Its central scintillation volume 
consists of 278 1 of ultra-pure PC (1,2,4-trimethylbenzene) doped with 1.5 g/1 of the fluor 
PPO (2,5-diphenyloxazole), contained in a spherical transparent 8.5 m diameter nylon Inner 
Vessel (IV). It is shielded by two buffer layers consisting of PC with the light quencher DMP 
(dimethylphthalate) . The surrounding 13.7 m diameter Stainless Steel Sphere (SSS) holds 
2212 inward-facing 8" photomultiplier tubes (PMTs) that detect scintillation light from the 
central region. The ID is surrounded by a powerful muon Outer Detector (OD) [15], composed 
of a high domed steel tank of 18 m diameter and 16.9 m height filled with 2 100 1 of ultra-pure 
water and instrumented with 208 PMTs that detect the Cherenkov emission from cosmic ray 
muons. 

3 The Cosmic Muon Flux 

The Gran Sasso underground laboratory consists of three experimental Halls, named A, B 
and C. These are about 80 m distant from each other and, in principle, could feature slightly 
different rock coverage. Borexino is located in Hall C, while the cosmic muon flux has been 
measured previously by LVD in Hall A [4] and by MACRO in Hall B [16]. These experiments 
reported a flux of (3.31 ± 0.03) • 10~ 4 ur V 1 and (3.22 ± 0.08) • 10 _4 m- 2 s _1 , respectively. In 
both cases the acceptance is strongly dependent on the muon incidence angle, contributing 
significantly to the systematic error of the final result. When comparing available results it 
should also be noted that the flux measured in different years can reflect differences in the 
mean air temperature (see table 1). 

This analysis is based on the first four years of Borexino data, taken between 16 th 
of May 2007 and 15 th of May 2011. The CNGS (CERN Neutrino to Gran Sasso [17, 18]) 
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Figure 1. Sketch of the Borexino detector. 



neutrino beam introduces muon events that are a background for this analysis. Consequently 
all events in coincidence with the beam spills are discarded (see [15] for details). The effective 
data set shows no prolonged or unevenly distributed downtime and accounts for a live time 
of 1063 days. Muon detection in Borexino is performed with both the ID and OD, however 
here we disregard events that generated a signal in the OD only. In this analysis the relevant 
cross-section for cosmic muons is therefore given by that of the SSS (146m 2 ), independent 
of the muon incoming angle. The corresponding total exposure of the data set is ~1.55T0 5 
m 2 -d and includes a sample of ~4.6T0 6 muons. 

The muons passing through the inner detector can be identified via three different 
methods. The first two are based on the detection of the Cherenkov light produced in the 
water: the light triggers the OD sub-system (MTF) or a cluster of hits is recognized within 
the time distribution of OD PMT hits (MCF). The third method (IDF) relies on the pulse 
shape identification of muon tracks among the point-like scintillation events detected by the 
ID. The detection efficiencies are 0.9925(2), 0.9928(2) and 0.9890(1) respectively. Details on 
the muon tagging methods and on how the efficiencies have been evaluated can be found in 
[15]- 

We have measured the muon rate through the ID using all strategies at our disposal and 
achieved consistent results. The average muon rate is (4310 ± 2 sta t ± 10 sys t) d _1 where the 
systematic error reflects the uncertainty in the efficiency and possible threshold effects. The 
rate corresponds to a cosmic muon flux of (3.41 ± 0.01) • 10~ 4 m _2 s _1 , taking into account 
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also the uncertainty in the SSS radius. This is the first measurement performed in Hall C 
and the first obtained with a spherical detector at LNGS. 



4 The Flux Modulation 



Air temperature increases during summer which lowers the average gas density. The less 
dense medium allows a longer mean free path of the mesons and increases the fraction of 
them that decay to produce muons before their first interaction. As only these muons are 
energetic enough to traverse the rock coverage of an underground site, a correlation between 
the muon flux observed underground and the air temperature is expected. We demonstrate 
such a correlation for the case of Borexino in section 8. Temperature fluctuations can have 
maxima and minima that occur at different dates in successive years and short term effects 
that are expected to perturb the ideal seasonal variation. Therefore a simple sinusoidal 
behavior is to be considered only a first order approximation. 

The muon flux measured day-by-day in Borexino is shown in figure 2 (upper panel) for 
the 1329 days for which valid data were available. A modulation is clearly visible. Fitting 
the distribution with the following function: 

I„ = ij + AI„ = ij + <5/ M cos {^{t - t )) (4.1) 

we obtain an average intensity If] = (3.414 ± 0.002 sta t) ■ 10 m~ 2 s , consistent with the 
flux reported in section 3, a period T = (366 ± 3) days, a modulation amplitude 51^ = 
(4.4 ± 0.2) • 10~ 6 m~ 2 s _1 , corresponding to (1.29 ± 0.07)% and a phase to = (179 ± 6) days, 
corresponding to a maximum on the 28 th of June; the Neyman's x 2 /NDF is 1558/1325. 

An alternative approach is to project and average the four years data set into one single 
year, as shown in figure 3. Fitting again with eq. 4.1 but with the period fixed to one year, 
we obtain consistent rate and amplitude. The phase is to = (179 ± 3) days. The x 2 /NDF of 
the fit is 442/362. 



5 The Atmospheric Model 

Deviations from the average muon flux that is measured underground, AJ jU (t) = — I®, 
can be related to variations from the average atmospheric temperature at a given altitude 
X, AT(X,t) = T(X,t) — T°(X) (from [6]). Considering every altitude layer, the net effect 
can be written as: 

poo 

AI M = / dXW(X)AT(X) (5.1) 
Jo 

where W(X) (see appendix A) reflects the altitude dependence of the production of mesons in 
the atmosphere and their decay into muons that can be observed underground. The integral 
extends over atmospheric depth from the altitude of muon production to the ground. 

The atmosphere can be described by many layers with a continuos distribution of tem- 
perature and pressure. A possible parametrization ([6] and with more details [19]) considers 
the atmosphere as an isothermal body with an effective temperature, T e g, obtained from a 
weighted average over atmospheric depth: 

f™dXT(X)W(X) Zn=o^X n T(X n )W(X n ) 
! eff = roo 7^7777^ - ^v . ~ 7Z77ZZ I I 5 - 2 ) 



J °° dXW(X) Zn=0^ n W(X n 
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Figure 2. Upper panel: cosmic muon signal measured by Borexino as a function of time. Lower panel: 
effective temperature, T e ff, computed using eq. 5.2 and averaging over the four daily measurements. 
Daily binning is used in both panels. The curves show the sinusoidal fit to the data (see text). 



xi<r 



10 

CM 

E 


0.355 


i] xn 


0.35 






c 




o 

3 


0.345 


5 






0.34 




0.335 




0.33 




0.325 



01/01 



02/03 



01/05 



01/07 



31/08 



31/10 



Figure 3. Cosmic muon flux: four years data set folded onto a one year period. Daily binning. The 
curve shows the sinusoidal fit to the data (see text). 



where the approximation may be done considering that the temperature is measured at 
discrete atmospheric levels, X n . 

Figure 4 shows the temperature in the atmosphere for the LNGS site and the weight 
function, W, as functions of the pressure levels. As can be seen, the higher layers of the 
atmosphere are given higher weights, as it is in these layers that most of the muons energetic 
enough to reach underground sites are produced. Muons produced at a lower altitude will 
be on average less energetic and a larger fraction of them will lie below threshold (E'thr)- 
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Figure 4. Average temperature (solid red line) [13] and normalized weight W[X) (black dashed line) 
as a function of pressure levels computed at the LNGS site. The right vertical axis shows the altitude 
corresponding to the pressure on the left vertical axis. 



We may also define the "effective temperature coefficient", ax, which quantifies the 
correlation effect that is discussed in section 8: 

a T = ^f dXW(X) (5.3) 
V Jo 



such that Eq. 5.1 may be written: 



~Jo = ( 5 - 4 ) 



6 Temperature Modulation 

The temperature data was obtained from the European Center for Medium-range Weather 
Forecasts (ECMWF)[13] which exploits different types of observations (e.g. surface, satellite, 
and upper air sounding) at many locations around the planet, and uses a global atmospheric 
model to interpolate to a particular location. In our case, the precise coordinates of the 
LNGS underground halls have been used: 13.5333° E, 42.4275° N. Atmospheric temperature 
is provided by the model at 37 discrete pressure levels in the [1-1000] hPa range (1 hPa = 
1.019 g/cm 2 ), four times a day at 00.00 h, 06.00 h, 12.00 h, and 18.00 h x . Based on this data 



1 The analysis in [3] and [4] used data from the air soundings performed by the Aeronautica Militare Italiana 
(AM) [20] near the military base of Pratica di Mare (12.44° E, 41.65° N), about 130km away from the lab. 
Aside to referring to a somewhat different location, that data set — probably the best available at the time of 
[3] — is significantly incomplete if compared to the one from ECMWF, both for number of measurements and 
for atmospheric depth coverage. We therefore used this data set only as a cross-check of the analysis based 
on the ECMWF data set, yielding consistent results. 
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set, T e g was calculated using eq. 5.2 four times a day. The four results were then averaged, 
and the variance of the four values was used to estimate the uncertainty in the mean. 

Figure 2 (lower panel) shows the daily values of T e fj for the four year period considered. 
A simple average gives T® s = 220.99 K, while the fit with a function analogous to eq. 4.1 
returns T e ° ff = (221.153 ± 0.007) K, amplitude (2.98 ± 0.01) K corresponding to 1.35 %, period 
(369.2 ± 0.2) days and phase (174.0 ± 0.4) days. These are very similar to the corresponding 
fit results of the muon flux data set, discussed in section 4. The x 2 /NDF = 46010/1457 
confirms that the sinusoidal behavior is only a first order approximation. Aside from small 
scale fluctuations, additional winter maxima can be observed which can be ascribed to the 
known meteorological phenomenon of the Sudden Stratospheric Warmings (SSW [21]) and 
whose effect is sometimes comparable in amplitude with the underlying seasonal modulation. 
In order to disentangle the seasonal dependence from sub-leading effects, the method of 
Lomb-Scargle has been used. 

7 Lomb-Scargle analysis 

Lomb-Scargle (LS) periodograms [22, 23] are a common method to analyze a binned data 
set for periodical modulations of the type 

N(t) =N -(1 + A- sm(2nt/T + <p)) . (7.1) 

Here, N(t) is the expected event rate at time t, while Nq represents the mean rate, A indicates 
the relative amplitude of the modulation, T describes its period and <p the phase relative to 
the start of the measurement. The power, P, for a particular modulation period, T, is defined 
as the weighted difference between the number of events, N(ti), in every data bin, i, and the 
weighted mean value, Nq, with cosine and sine functions that oscillate with an investigated 
period, T: 

p= 1 f lEtiMNjtd-N^cosTi} 2 | [^ =1 w l (N(t l )-N )sm n } 2 

where n is the number of bins and t{ the time at which the data corresponding to bin i was 
acquired. The weights Wi = a^ 2 / (cr^ 2 ) are the inverse squares of the statistical uncertainties 
of individual bins, N(ti), divided by their average value. Accordingly, a 2 is the weighted 
variance of the data. The phase of the sine and cosine weights, Tj, is defined as [24] 

with tr = l^J Zl^ffi) \. ( , 3) 

As the quadratic sums of both cosine and sine are used in eq. (7.2), the result is indepen- 
dent of the modulation phase as long as the modulation period is shorter than the overall 
measurement time. 

Figure 5 shows an LS periodogram of the four-year muon data acquired in Borexino. 
The LS power, P, of a given modulation is primarily a function of its amplitude, A. Statistical 
fluctuations of the bin content will alter both the maximum, P, generated by white noise and 
the exact value observed for the actual modulation. To assess the significance of a modulation 
discovery, it is therefore necessary to know the statistical fluctuations of both the white noise 
level and the signal height. 
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Figure 5. LS periodogram of the muon data. 
The dashed line indicates the detection threshold 
(3a). 
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Figure 6. LS periodogram of the temperature 
data. The dashed line indicates the detection 
threshold (3cr). 



Using the known detector live time distribution and muon rate, we have performed 
realistic Monte Carlo (MC) simulations of muon sample time distributions. In a first step, 10 4 
MC data samples with a constant muon rate have been simulated. From the corresponding 
P spectra we can estimate the probability that an apparent periodic modulation appears as 
a result of white noise. For a given period T, the 3a detection threshold P t hr(T) is set so 
that 99.7% of all white noise samples feature lower values of P. This threshold is indicated 
by the dashed line in figure 5. Over a broad range of periods, we obtain values of P t h r ~ 7 

(3ff). 

Based on a further set of MC data samples that contain periodic modulations, given 
values of P can be matched to modulation amplitudes and vice versa: Pthr corresponds to a 
relative amplitude of ~0.3 %, the exact value depending on statistical fluctuations. Compared 
to threshold, the P peak corresponding to the annual modulation is highly significant. The 
maximum of P = 140 is reached for a period of 364 days. From MC simulation, one can 
associate this LS power to an amplitude of (1.20 ± 0.05) %. 

In figure 5, several secondary peaks are visible: two peaks at T = 0.05 years and T = 0.5 
years are just above detection threshold, while the largest secondary peak is at T = 1.7 years 
featuring P ~ 20, corresponding to an amplitude of about 0.4%. The appearance of these 
additional peaks indicates that there are deviations of the modulation pattern from a simple 
year-long sinusoidal. Dividing the data sample by the main modulation allows to investigate 
whether the side-peaks are reflections of the annual modulations at harmonic frequencies. 
This procedure increases the significance of the peak at T = 0.5 years to P = 15, while both 
the peaks at T = 0.05 years and T = 1.7 years are pushed below detection threshold. To 
further investigate this semi-annual modulation, a second sinusoid was added to the direct 
fit to the data set (eq. 4.1). The period of the second sinusoid fits to T = (179 ± 2) days, 
and the amplitude to A = (0.37 ± 0.07) %. The fit result for the main (1-year) oscillation 
remains unaffected within uncertainties. We conclude that the modulation found in the 
data is best described by the superposition of two sinusoidal terms with the semi-annual 
component slightly distorting the rising and falling flanks of the main annual component. 

The LS analysis can also be applied to the temperature data. The resulting periodogram 
of figure 6 features again a prominent peak at T = 368 days. The observed value of P = 367 
is compatible with the modulation amplitude of 1.35% (section 6) with 90% confidence 
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level. The observation threshold (3a) for the temperature data set is P t h r (ly) = 58; none 
of the secondary peaks is significant at 3a level. The peak corresponding to the semi-annual 
sub-modulation identified for the muons has a ~2.2cr significance in the temperature data. 



8 Correlation 



Figure 2 shows the correlation between fluctuations in the atmospheric temperature and the 
cosmic muon flux. To quantify such a correlation we plot AI^/I® vs AT c ff/T^ for each 
day in figure 7. Only days with duty cycle > 50% are included for a total of 1165 days. 
The correlation coefficient (R-value) between these two distributions is 0.60 indicating a 
positive correlation. We now want to determine the effective temperature coefficient (eq. 5.3). 
We perform a linear regression accounting for error bars on both axes using a numerical 
minimization method. As a result we obtain olt = 0.93±0.04 sta t with x 2 /NDF = 1144/1164. 
This result is consistent and features smaller errors when compared to aj 1 = 0.91 ± 0.07, the 
previous measurement by MACRO at Gran Sasso [25]. 

We have performed the following tests to check for systematic uncertainties: 

• If] and T e ° ff have been computed in different ways: averaging 1^ and T e g values over the 
available data set; from the fit to the four year data set with the sinusoidal functions as 
in eq. 4.1 and figure 2; from a fit of the same data set with a constant function. In addi- 
tion T® s has been computed including or excluding the days for which no corresponding 
muon flux was available. 

• The analysis has been performed on a moving two-year sub-period checking the stability 
of the result. 
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Figure 8. Measured values for the effective temperature coefficient, ar, at varying site depths. The 
results from this analysis (in blue) as well as those from different experiments are presented. The red 
line is the value predicted including muon production by pions and kaons. The dashed lines account 
for one production mechanism only. See [6] and refs therein for details. 



• We have varied the requirement of including only days with duty cycle > 50% in the 
range [> 20%, > 80%]. 

• We have considered E^hr = 1.833 TeV from [19] and £thr = 1.3 TeV as in [3, 4] for the 
computation of T e g (see also appendix A). 

• We ran the analysis substituting the air temperature data set from ECMWF with that 
from Aer. Mil. Italiana (see footnote in section 6) used in [3, 4]. 

• We recomputed T e g using weights that account for muon production only from pion 
decay, i.e. neglecting the kaon contribution as done in [3, 4] (see appendix A). 

In all cases we found our result to be stable, so we believe that the systematic uncertainty is 
small compared to the statistical error from the fit. 

Figure 8 shows the measured values for ay, with the details summarized in table 1. 
This value has been predicted for different site depths following [6]. As shown in figure 8, 
otT asymptotically approaches unity with increasing site depth. This is because the air- 
density-independent contribution to the muon signal originating from mesons which have 
interacted before decaying is progressively left below threshold. At LNGS ay is expected to 
be 0.92 ± 0.02 (considering muon production from both pions and kaons) in good agreement 



Site [km w.e.l 

L J 


S. Pole(1.6) 


Soudan(2.1) 


LNGS-A(3.8) 


LNGS-B(3.8) 


LNGS-C(3.8) 


Experiment 


ICECUBE[5] 


MINOS [6] 


LVD [4] 


MACRO [25] 


Borexino 


£ t hr[GeV] 


466 


730 


1833 


1833 


1833 


Rate [10 3 /i/d] 


100000 


40 


8 


9 


4 


Meas. time [y] 


2007-ll[4y] 


2003-08 [5y] 


2001-08[8y] 


1991-97[7y] 


2007-ll[4y] 


Accept. [cos(#)] 


any 


[0.05-0.92] 


>0.5 


>0.3 


any 


UW-^m-Zs- 1 ] 






3.31 ±0.03 


3.22 ±0.08 


3.41 ±0.01 


Modul. Ampl. 






1.5% 




1.3% 


Period (days) 






367 ± 15 




366 ±3 


Phase (days) 






185 ± 15 
(Jul 4 th ) 




179 ±6 
(Jun 28 th ) 


Binning 


daily 


daily 


daily 


monthly 


daily 


Air Data 


NASA-AIRS 


ECMWF 


Aer.Mil. 


Aer.Mil. 


ECMWF 


T e ff model 


7T+K 


7T+K 


7r-only 


7r-only 


7T±K 


Correlation 


n.p. 


0.90 


0.53 


0.91 


0.62 


corneas.) 


0.860 ±0.010 


0.873 ± 0.009 




0.91 ±0.07 


0.93 ± 0.04 


ar(pred.) 


~ 0.83 


0.864 ±0.024 


0.92 ±0.02 


0.92 ±0.02 


0.92 ±0.02 



Table 1. Comparison of the different analyses of the muon seasonal modulation and correlation with 
temperature by some existing underground experiments. 



with the result from this analysis. The systematic uncertainty in this prediction was found 
by modifying the input parameters for the computation according to their uncertainties and 
recalculating. 

With a longer exposure we expect to measure ax with better precision, opening the 
way to an indirect determination of the K/ir ratio in the interaction of primary cosmic rays 
in the atmosphere with the method detailed in [6, 19] and probing a complementary energy 
region compared with existing accelerator experiments. 

9 Conclusions 

Borexino has reached four years of continuous data taking at LNGS under a rock coverage of 
3800 mw.e. Due to the spherical geometry of the detector, we have measured the underground 
cosmic muon flux with reduced systematics: (3.41±0.01)T0 _4 m _2 s _1 . We also have observed 
a seasonal modulation of the flux and measured the amplitude to be (1.29 ± 0.07) % and the 
phase to be (179 ± 6) days corresponding to a maximum on the 28 th of June. To invesitgate 
the correlation between air temperature variations and changes in the muon flux, we have 
obtained air temperature data from global atmospheric models for the same four years period. 
We showed that the annual modulation is also present in the effective temperature data, 
with oscillation parameters compatible with those of the muon modulation. We then showed 
that the two data sets are positively correlated (R = 0.60) and we measured the effective 
temperature coefficient to be qt = 0.93 ± 0.04. This result is compatible with theoretical 
expectations and is an improvement in precision from previous measurements at Gran Sasso. 
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A Effective temperature weight function 

The weight W{X) used in eq. 5.2 can be written as the sum W n + Wk, representing the 
contribution of pions and kaons to the overall variation in muon intensity: 

W*' K (X) ~ (1 ~ X/K'^ tK ) 2 e~ x I^^Al^ K 

K > 1 +( 1 + l)Bl K K(X)((E thl co S e)/e 7r , K y { ^ 

where 

(1 - X/AL K f 

K{X) = V vl ' W ' K) (A.2) 

(l-e- X ^)K <K /X 

The parameters A\ K include the amount of inclusive meson production in the forward frag- 
mentation region, the masses of mesons and muons, and the muon spectral index 7; the input 
values are A\ = 1 and A l K = 0.38-r K / n , where r K / w is the K/ir number ratio. The parameters 
B\ k reflect the relative atmospheric attenuation of mesons; the threshold energy, -Ethr, is the 
energy required for a muon to survive to a particular underground depth and 9 is the angle 
between the muon and the vertical directions; the attenuation lengths for the cosmic ray 
primaries, pions and kaons are A^r, A w and Ajc respectively with l/A' w K = 1/Ajv — I/A^k- 
The meson critical energy, e^ Ki is the meson energy for which decay and interaction have an 
equal probability. The value of (Ethr cos 9) used here is the median of the distribution. The 
values for these parameters can be found in table 1 of [6], with the exception of (i? t hrcos0) 
which is site dependent and is found by MC simulations. At LNGS {E t \ iV cos9) = 1-833 TeV 
according to [19]. The dependence of W(X) on E t hr is however moderate. 
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